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Abstract 



We present results for the light quark masses obtained from a lattice QCD simulation 
with Nf — 2 degenerate Wilson dynamical quark flavours. The sea quark masses of our 
lattice, of spacing a ~ 0.06 fm, are relatively heavy, i.e., they cover the range corresponding 
to 0.60 < MpjMy < 0.75. After implementing the non-perturbative RI-MOM method to 
renormalise quark masses, we obtain m^f(2GeV) = 4.3 ± 0.4T.Q 1 MeV, and m^ s (2GeV) = 
lOliS^Q 5 MeV, which are about 15% larger than they would be if renormalised perturbatively. 
In addition, we show that the above results are compatible with those obtained in a quenched 
simulation with a similar lattice. 



1 Introduction 



An accurate determination of quark masses is highly important for both theoretical studies of 
flavour physics and particle physics phenomenology. Within the quenched approximation (Nj = 0), 
lattice QCD calculations reached a few percent accuracy, so that the error due to the use of 
quenched approximation became the main source of uncertainty. To examine the effects of the 
inclusion of the sea quarks several lattice results with either Nf = 2 or Nf = 3 dynamical quarks 
have recently been reported [1-9]. Many of them seem to indicate that the unquenching lowers the 
physical values of light quark masses. For example, the strange quark mass, which in the quenched 
approximation is about m^ s (2GeV) ~ 100 MeV, becomes m^ s (2GeV) ~ 85 4 90 MeV when 
Nf = 2 dynamical quarks are included [2, 3], and even smaller with Nf = 3, namely m^ s (2 GeV) ~ 
75 4-80 MeV [6,7]. 

The accuracy of the results obtained from simulations with Nf ^ 0, however, is still limited by 
several systematic uncertainties. In particular, precision quenched studies of light quark masses 
evidenced the importance of computing the mass renormalisation constants non-perturbatively. 
In almost all the studies with dynamical fermions performed so far, however, renormalisation has 
been made by means of one-loop boosted perturbation theory (BPT). Only very recently, non- 
perturbative renormalisation (NPR) has been implemented in a simulation with Nf = 2 by the 
QCDSF collaboration. The resulting strange quark mass value is m^ s (2GeV) = 119±9 MeV [8], 
thus larger by about 20% than the one previously reported in ref. [3], where the same lattice action 
and the same vector definition of the bare quark mass have been used, but with the renormalisation 
constant estimated perturbatively. Recently, also the Alpha collaboration reported their value for 
the NPR-ed strange quark mass, m^ s (2GeV) = 97 ± 22 MeV [9]. It appears that the results of 
both QCDSF and Alpha are consistent with their previous quenched estimates [10, 11]. 

In this paper we present our results for the light quark masses obtained on a 24 3 x 48 lattice 
with Nf = 2 dynamical mass-degenerate quarks, in which we implement the RI-MOM NPR 
method of refs. [12, 13]. In the numerical simulation we used the Hybrid Monte Carlo algorithm 
(HMC) [14, 15], and chose to work with the Wilson plaquette gauge action and the Wilson quark 
action at (3 = 5.8. Though the fermionic action is not improved at 0(a), this value of the gauge 
coupling corresponds to a rather small value of the lattice spacing, namely a ~ 0.06 fm. The 
spatial extension of our lattice is L ~ 1.5 fm. Gauge configurations have been generated with four 
different values of the sea quark mass, for which the ratio of the pseudoscalar over vector meson 
masses lies in the range Mp/My ~ 0.60 4 0.75. Our final results for the average up/down and 
strange quark masses are 

m5|(2GeV)=4.3(4)(+i- 1 ) MeV / N f = 2 

m^(2GeV) = 101 (8) (If) MeV \ RI-MOM 

where the first error is statistical and the second systematic, the latter dominated by discretisation 
lattice artifacts. Our central values refer to the quark mass definition based on the use of the axial 
Ward identity. The above results are larger than the ones we obtain by using one-loop BPT to 
compute the mass renormalisation, namely 

m]f(2GeV) = 3.7(3)(^- 3 )MeV / N f = 2 

mp(2GeV) = 88(7) (t 30 ) MeV ^ BPT 

The latter values agree with the Nf = 2 results reported in refs. [2, 3] in which BPT renormalisation 
has been used, while the values in eq. (1) agree with refs. [8, 9] in which NPR has been implemented. 
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In performing such a comparison it should be also kept in mind, however, that the results of 
refs. [2, 3, 8, 9] are affected by discretization errors of different size, having been produced by using 
C(a)-improved actions but at larger values of lattice spacings. In addition, the results of refs. [2, 8] 
have been obtained after performing an extrapolation to the continuum limit. The values given in 
eq. (1) can be also compared with our quenched estimate, 



m^ b (2 GeV) = 4.6(2)(+^ & ) MeV / N f = 

1^(2 GeV) = 106(2) 



mf s (2 GeV) = 106(2) (+J 2 ) MeV V RI-MOM 



(3) 



obtained on a lattice similar in size and resolution, and by employing the same Wilson quark 
action. We therefore conclude that, within our statistical and systematic accuracy, and with the 
sea quark masses used in our simulations, we do not observe any significant effect due to the 
presence of dynamical quarks. We plan to pursue this study by exploring lighter sea quarks. 

The remainder of this paper is organised as follows: in sec. 2 we discuss the details of our 
numerical simulations; in sec. 3 we present the results for the meson masses and the bare quark 
masses, directly accessed from the simulation; in sec. 4 we provide the mass renormalisation factors, 
while in sec. 5 we determine the physical quark mass values; we briefly conclude in sec. 6. 

2 Simulation details 

In our numerical simulation we choose to work with the Wilson action 

S w = S g + S q , (4) 
where S g is the standard Wilson plaquette gauge action, 

s * = § E Tr u * = § E Tr p^ u *+»M + ^ u U > ( 5 ) 

p X,/J,U 

and Nf = 2 flavours of mass-degenerate dynamical quarks are included by using the Wilson action 



S q — QxD xy q y — Qx 

x,y x,y 



S xy ~ - iJUx^Sx+hy + (1+7^)^4,^,2/+^] 



% , (6) 



with k being the usual hopping parameter. The simulation parameters used in this work are 
summarised in table 1. To generate the ensembles of gauge field configurations we employ the 
HMC algorithm [14, 15] and implement the LL-SSOR preconditioning [16]. In solving the molecular 
dynamics equations we use the leapfrog integration scheme. Each trajectory of unit length (r = 1) 
is divided into N step = 250 time-steps. The resulting St = T/N stcp = 4 x 10~ 3 appears to be 
small enough to avoid large energy differences along a trajectory, thus increasing the acceptance 
probability. We find an acceptance probability of about 85%. To check the reversibility along 
the trajectory we verify that the quantity — 1|| remains unchanged at the relative level of 

10~ n . For the quark matrix inversion algorithm we used the BiConjugate Gradient Stabilized 
(BiCGStab) method [17]. The inversion rounding errors are of the order of 10~ 8 . 

To select sufficiently decorrelated configurations we studied the autocorrelation as a function 
of the number of trajectories t by which two configurations are separated. For a given observable 
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(3 = 5.8 


f3 = 5.6 


K s 


n 1 rnr fi i ron A 1 £ /I A A 1 C A 1 

U.looo U.lOoO U.104U U.1041 


A 1 C£!A A 1 A 1 CCQA 

U.iOOU U.iOiO U.iOSU 


L 3 xT 


24 3 x 48 


16 3 x 48 


minutes/traj. 
Tmc 

N C onf 
(P) 


10 12 16 21 
2250 2250 2250 2250 
50 50 50 50 
0.59268(3) 0.59312(2) 0.59336(3) 0.59341(3) 


5.3 13 28 
2250 2250 2250 
50 50 50 
0.56990(6) 0.57248(6) 0.57261(6) 


L 3 x T 


16 3 x 48 




minutes/traj. 
Tmc 

Nconf 
(P) 


3.3 4.4 4.4 5.6 
4500 4500 4500 4500 
100 100 100 100 
0.59283(3) 0.59314(3) 0.59340(3) 0.59345(3) 



Table 1: Run parameters for the simulations with Nf = 2. From the entire ensemble of generated HMC 
trajectories of unit length, Tmc, w & select N con f for the data analysis. The stopping condition in the 
quark matrix inversion (D xy G y ~ B x ) is \\DG — B\\ < 10~ 15 . (P) is the value of the average plaquette. 
"Minutes/traj. " refers to the time required to produce one trajectory on the APEmillc machine (128 GFlop). 



A, the autocorrelation is evaluated as 

-. T M c—t 

c ' 4 ( ( > = frss ■ < 7) 

Tmc „i 

where (A) is the mean value of A computed on the total number of trajectories Tmc and A s 
denotes the value measured along the trajectory s. After examining the autocorrelation of the 
plaquette and of the pseudoscalar two-point correlation functions, both shown in fig. 1, we made 
the conservative choice to select configurations separated by 45 trajectories. 1 The sea quark 
masses used in our simulations correspond to ratios of pseudoscalar over vector meson masses 
covering the range 0.60 < Mp/My < 0.75. In our data analysis we distinguish the sea (k s ) from 
the valence (n v ) hopping parameters, where, for each value of the sea quark mass, k v can assume 
any of the values listed in table 1. The main results of this paper refer to the run made at f3 = 5.8 
and L = 24 (i.e., L ~ 1.5 fm), whereas those obtained on the lattice with L = 16 and (3 either 
5.8 (L ~ 1 fm) or 5.6 (L ~ 1.3 fm) are used to investigate finite volume and discretisation effects 
respectively. 

Finally, in order to assess the effects of quenching, we also made a quenched run (250 config- 
urations) at (3 = 6.2 and V = 24 3 x 56 (L ~ 1.6 fm). To reproduce the values of the Nf = 2 
pseudoscalar meson masses, we chose in the quenched run k v £ {0.1514,0.1517,0.1519,0.1524}. 

An important feature of the data analysis in the partially quenched case is the appropriate 
treatment of the statistical errors. While the results of computations of an observable O on 
configurations with different values of k s are statistically independent, those computed on the 
gauge configurations with the same ft, but different k v are not. To account for these correlation 



1 The autocorrelation of the pseudoscalar two-point function is made for the time separation between the sources 
fixed to 14. 
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PLAQUETTE AUTOCORRELATION PSEUDOSCALAR AUTOCORRELATION 




25 50 75 100 125 150 25 50 75 100 125 150 

t t 

Figure 1: Autocorrelations of the plaquette C p (t) (left), and of the pseudoscalar correlation function C p (t) 
(right), as a function of the number of trajectories t used to separate consecutive configurations. For C PS (t) 
the fluctuations around zero start at about t = 30, indicating that conhgurations separated by less than 30 
trajectories are correlated. The above plots refer to the run at (3 — 5.8 (V — 24 3 x 48, and n s = 0.1538). 



effects we use a combination of the bootstrap and jackknife statistical methods. At fixed k s , our 
data are clustered into Njk subsets and the statistical errors on O are estimated by using the 
jackknife method. The resulting Njk values at a given k s , O^, are then considered as a set of 
uniformly distributed numbers. In a next step one generates iVj, S> Njk bootstrap events by taking, 
for each n s , one of the uniformly distributed Ofy. The variance of O is then estimated by 

o- 2 = (N jk - 1) x f(0 2 ) Nb - (0) 2 Nl ) , (8) 

where {--^Nt stands for the average calculated over the Nb bootstrap events, while the factor 
(Njk — 1) takes into account the jackknife correlation at fixed k s . For a comparison, we also 
considered the procedure adopted in ref. [2], consisting in applying the jackknife method at one 
specific value of k s , while keeping the data at the remaining three k s fixed at their central values. 
As expected, after adding in quadrature the four variances calculated in this way, we find excellent 
agreement with the estimates obtained by using eq. (8). 

3 Masses of light mesons and quarks 

In this section we present the results for the light pseudoscalar and vector meson masses, as well as 
the bare quark masses, obtained from our main data-set (/3 = 5.8, L = 24). These quantities are 
extracted from the standard study of two-point correlation functions, Cjj(t), in which we consider 
both degenerate and non-degenerate valence quarks, with valence quarks either equal or different 
from the sea quarks. At large Euclidean time separations, f > 0, Cjj(t) is dominated by the 
lightest hadronic state and we fit it to the form: 

Cjj(t) = J2(0\J(S,t)JH0)\0) m |^( e -^ + r?e -^M) , (9 ) 

X 
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K Vl 


K V2 


M P 


My 


m^(a) 


0.1535 


0.1535 


0.1535 


0.262(4) 


0.348(8) 


0.0333(6) 




0.1535 


0.1538 


0.252(5) 


0.341(8) 


0.0306(6) 




0.1535 


0.1540 


0.245(5) 


0.336(9) 


0.0289(6) 




0.1535 


0.1541 


0.241(5) 


0.334(9) 


0.0280(6) 




0.1538 


0.1538 


0.241(5) 


0.333(9) 


0.0280(6) 




0.1538 


0.1540 


0.233(5) 


0.329(9) 


0.0263(6) 




0.1538 


0.1541 


0.229(5) 


0.326(9) 


0.0254(6) 




0.1540 


0.1540 


0.226(5) 


0.324(9) 


0.0246(6) 




0.1540 


0.1541 


0.222(5) 


0.322(10) 


0.0237(6) 




0.1541 


0.1541 


0.218(5) 


0.319(10) 


0.0228(6) 


0.1538 


0.1535 


0.1535 


0.258(3) 


0.335(4) 


0.0314(2) 




0.1535 


0.1538 


0.247(3) 


0.328(4) 


0.0287(2) 




0.1535 


0.1540 


0.240(3) 


0.323(5) 


0.0269(2) 




0.1535 


0.1541 


0.236(4) 


0.321(5) 


0.0261(2) 




0.1538 


0.1538 


0.236(4) 


0.321(5) 


0.0260(2) 




0.1538 


0.1540 


0.228(4) 


0.316(6) 


0.0243(2) 




0.1538 


0.1541 


0.224(4) 


0.314(6) 


0.0234(2) 




0.1540 


0.1540 


0.220(4) 


0.311(6) 


0.0226(2) 




0.1540 


0.1541 


0.216(5) 


0.309(7) 


0.0217(2) 




0.1541 


0.1541 


0.212(5) 


0.306(7) 


0.0208(2) 


0.1540 


0.1535 


0.1535 


0.258(2) 


0.345(5) 


0.0303(3) 




0.1535 


0.1538 


0.247(2) 


0.338(6) 


0.0277(3) 




0.1535 


0.1540 


0.240(3) 


0.333(7) 


0.0259(3) 




0.1535 


0.1541 


0.236(3) 


0.331(7) 


0.0250(3) 




0.1538 


0.1538 


0.236(3) 


0.331(7) 


0.0250(3) 




0.1538 


0.1540 


0.228(3) 


0.326(8) 


0.0233(3) 




0.1538 


0.1541 


0.225(3) 


0.324(8) 


0.0224(3) 




0.1540 


0.1540 


0.221(3) 


0.321(9) 


0.0215(4) 




0.1540 


0.1541 


0.217(3) 


0.319(9) 


0.0207(4) 




0.1541 


0.1541 


0.212(3) 


0.316(10) 


0.0198(4) 


0.1541 


0.1535 


0.1535 


0.234(6) 


0.318(10) 


0.0286(3) 




0.1535 


0.1538 


0.222(6) 


0.311(11) 


0.0259(3) 




0.1535 


0.1540 


0.213(6) 


0.307(11) 


0.0242(3) 




0.1535 


0.1541 


0.209(6) 


0.304(11) 


0.0233(3) 




0.1538 


0.1538 


0.209(6) 


0.304(11) 


0.0233(3) 




0.1538 


0.1540 


0.200(7) 


0.300(12) 


0.0215(3) 




0.1538 


0.1541 


0.196(7) 


0.298(12) 


0.0206(3) 




0.1540 


0.1540 


0.191(7) 


0.295(12) 


0.0198(3) 




0.1540 


0.1541 


0.187(7) 


0.293(13) 


0.0189(3) 




0.1541 


0.1541 


0.182(7) 


0.291(13) 


0.0180(4) 



Table 2: Pseudoscalar and vector meson masses, and the bare quark masses m^ l (a) — |[m„i(a) + 

m V 2{a)] AWI obtained by using the axial Ward identity method (cf. eq. (10)). Results, in lattice units, refer 
to the simulation with /3 — 5.8 and V — 24 3 x 48. Those obtained at (3 = 5.6 are tabulated in appendix 1. 
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where rj is the time reversal (t <-> T — t) symmetry factor. 2 The pseudoscalar and vector me- 
son masses are extracted from Cpp(t) and Cyv(t) respectively, by fitting in t £ [14,23]. The 
corresponding results are listed in table 2. These results are obtained by using local source opera- 
tors. We checked, however, that the Jacobi smearing procedure [18] leaves the masses unchanged, 
although the lowest hadronic state becomes isolated earlier. 

In the same table 2 we also give the results for the quark masses obtained by using the axial 
Ward identity (AWI) definition. More specifically, with Aq(x) = <foi7o75(Zu2) and P = q v i75Qv2, 

(Y,d A o (x)P\0)} 

x <>0 r , s . , ^-|AWI , 1fV . 

— — > \m vl (a) + m v2 {a)\ , (10) 

(^P(x)pt(o)) 

X 

where we fit in the same interval as before, i.e., t £ [14,23]. Alternatively, one can use the vector 
Ward identity (VWI) and relate the bare quark mass to the Wilson hopping parameter as: 

with k v being either k v \ or k v2 (cf. table 2), and K cr (K s ) is the critical hopping parameter for 
a given sea quark mass. Notice that both the AWI and VWI definitions of the quark mass do 
depend on the values of the sea quark masses. In the axial case, this dependence is provided by the 
ratio of two-point correlation functions of eq. (10). In the VWI definition, the dependence on the 
sea quark mass is given by the value of the critical hopping parameter K cr {n s ) computed at fixed 
k s . We also notice that the sea quark dependence of the critical hopping parameter is numerically 
equivalent to the observation made in ref. [8] where such a dependence is expressed in terms of the 
renormalisation constants of the singlet and non-singlet scalar density. 



3.1 Lattice spacing and /t cr (/t s ) 

In order to fix the value of the lattice spacing and K cr (K s ), we inspect the dependence of the 
pseudoscalar and vector meson masses on the valence and sea quark masses. The left panel of 
fig. 2 shows the dependence of the vector meson masses on the squared pseudoscalar ones, which 
suggests a simple linear fit to the form 

M v (k s , k v1 ,k V2 ) = P 1 + P 2 x M p (k s , k v1 ,k V2 ) + P 3 x Mp(K s , k s , k s ) , (12) 

in an obvious notation, and with capital letters used to denote that the meson masses are given 
in lattice units. Such a fit yields P3 = 0.1(3), i.e., the sea quark dependence is very weak. We 
can therefore neglect it and set P3 = 0. In this way we obtain P\ = 0.25(1) and P 2 = 1-4(2). We 
then determine the lattice spacing from this fit and by using the method of the "physical lattice 
planes" [19]: at the intersection of our data with the physical value for the ratio mx/mx* = 0.554, 
we impose that Mk* obtained on the lattice is Mk* = am p ^* s thus obtaining 

a- 1 ^ =3.2(1) GeV. (13) 

We also estimate the value of the lattice spacing by studying the static quark potential. At each 
value of the sea quark mass we extract the following values of the Sommer parameter Rq = r^/a [20] : 

(Rq) Ks ={7.52(17) ai5 35 , 7.70(9)o.i538 , 7.78(19) .i 5 40 , 8.10(20) .i 5 4i}. (14) 

2 rj — +1 for Cpp(t), Cvv(t), or Caa(£), whereas for CVip(t), rj — —1. In this notation P — qjsq, V = Vi — qyiq, 
and A = A Q = 97075?- 
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Figure 2: Data-points that illustrate the functional dependencies discussed in eqs. (12) and (16) 



Details on this extraction can be found in appendix 2. Assuming the observed dependence of R$ 
on the sea quark mass to be physical, as the lattice spacing a is supposed to be independent on 
the sea quark mass, we linearly extrapolate the above results to the chiral limit [Rq = a + (3 x 
Mp(n s , k s , k s )], ending up with R = 8.6(4). After setting r = 0.5 fm, we find 



a* 1 = 3.4(2) GeV 



(15) 



in good agreement with the value given in eq. (13). 

Finally, for my wi in eq. (11) we also need k ct .(k s ). Inspired by partially quenched ChPT [21], 
but neglecting the chiral logs since our dynamical quarks are not light enough, we fit our data to 



M£(k s ,k v1 ,k V2 ) = Q 1 m v v wl [l + Q 2 m v s wl +Q 3 <H, 



(16) 



and obtain Q2 = —0.2(25) and Q3 = 0.4(14). In other words our data are not sensitive to the 
quadratic corrections proportional to Q2 and Q3, as it can be seen from the right plot in fig. 2. 
Very similar feature is observed if instead of 

we set Q2 = Q3 = 0. From the resulting fit to eq. (16), we get Qi = 1.70(3) and 



1 i!f we use rrivY 1 \ an d therefore in the following 



(«cr) Ks ={0.15544(7)o.i535, 0.15537(6)o.i538 , 0.155 37(5)o.i540, 0.15503(8)o.i54i } 



(17) 



4 Quark mass renormalisation 

Before discussing the strategy to identify the physical strange and the average up/down quark 
masses from the results obtained by using eqs. (10,11), we need to determine the corresponding mul- 
tiplicative mass renormalisation constants, Z^ WI (a/j,) = Z^/Zp(a/j,), and Z^ VI (a\i) = \ jZs{a\x). 
For completeness, we also present in this section the results for Zy, Z? and the quark field RC Z q . 

We use the RI-MOM method [12] which we explained in great detail in a previous publica- 
tion [13]. We adopt the same notation as in ref. [13] and compute the renormalisation group 
invariant combinations 

Z o (afi ) = C {a^) Zg GI = C o (afi ) (Z (an)/C (an)) , (18) 
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for the scale dependent bilinear operators, O = S,P,T. The evolution functions Co(a/x), which 
are known in the RI-MOM scheme at the N 3 LO for Z s and Z P [22] and at the N 2 LO for Z T [23], 
explicitly cancel the scale dependence of the RCs, Zo(afi), at the considered order in continuum 
perturbation theory. 

In fig. 3 we plot Zy,Ai as well as Zs,p,T(ano), where in eq. (18) we choose a/io = 1. All of them 
are supposed to be flat in (a/i) 2 , where /i is the initial renormalisation scale introduced in eq. (18), 
modulo discretisation effects. In fact, however, we see that these effects are pronounced in the 
case of ZsiflHo), which we correct for by linearly extrapolating from 1 < (a/i) 2 < 2, to (a/x) 2 = 0. 
Furthermore, for the RI-MOM scheme to be mass independent we extrapolate the mild sea quark 

0.85 



0.75 



0.65 




0.55 











♦ z/n„= 


l/a) 




1/a) 




l/a) 



0.0 



0.5 
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(a\if 



1.5 



2.0 



Figure 3: Renormalisation constants in the RI-MOM scheme at the scale a/xo = b computed according to 
eq. (18), as a function of the initial renormalisation scale {an) 2 . Discretisation effects cx (a/i) 2 are eliminated 
by extrapolating the data from 1 < (a/i) 2 < 2, to (a/i) 2 = 0. Illustration is provided for the data with 
k s = 0.1538. 



dependence linearly to the chiral limit. The results are reported in table 3. By confronting the 
results for the RCs obtained at /3 = 5.8 and L = 24 with those obtained at the same but with 
L = 16, we find that they differ by less than 5%. Given the smallness of the physical lattice size at 
L = 16 (~ 1 fm), we expect finite volume effects on our estimate of the RCs to be negligible with 
respect to discretization errors, which are evaluated to be of the order of 10% (see subsec. 5.1). 
This is expected, since the RCs encode the short distance physics while the volume finiteness affects 
the long distance physics. For an easier comparison of the non-perturbatively and perturbatively 
estimated RCs, in table 3 we also give the RCs evaluated by means of 1-loop BPT [24]. These are 
obtained by replacing in the one-loop expressions of ref. [25] the bare lattice coupling by a boosted 
one, (?o ~~ > ff 2 - We distinguish the so-called "naive" boosting 
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2 - al/{P) , 



(19) 



with (P) being the average plaquette, and the so called "tadpole improved" coupling, which is 
defined by inverting the perturbative series of the logarithm of the plaquette, namely 3 



ln(P) 



-^ 2 (3.40/a) 



1 - (1.1905 - 0.2266^/) 



5 2 (3-40/a) 



47T 



(20) 



! For a more extensive discussion about various forms of boosting, see ref. [26]. 
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P(N f ) 




z q 


Z v 


Z A 


Z s 


z P 




5.8 (N f = 2) 


0.1535 


0.81(1) 


0.70(1) 


0.80(2) 


0.74(1) 


0.54(1) 


0.78(2) 




0.1538 


0.80(1) 


0.70(1) 


0.79(2) 


0.73(2) 


0.55(1) 


0.77(2) 




0.1540 


0.78(1) 


0.67(1) 


0.77(2) 


0.72(1) 


0.52(1) 


0.75(2) 




0.1541 


0.79(1) 


0.69(1) 


0.78(2) 


0.71(2) 


0.52(1) 


0.77(2) 




K cr ( K cr ) 


0.76(1) 


0.66(2) 


0.76(1) 


0.67(3) 


0.50(2) 


0.74(4) 




N-BPT 


0.75 


0.70 


0.77 


0.75 


0.61 


0.75 




TI-BPT 


0.72 


0.65 


0.73 


0.71 


0.55 


0.71 


5.6 (Nf = 2) 


0.1560 


0.78(2) 


0.63(4) 


0.78(3) 


0.67(2) 


0.46(1) 


0.76(5) 




0.1575 


0.78(2) 


0.64(4) 


0.77(2) 


0.65(2) 


0.46(1) 


0.77(3) 




0.1580 


0.78(1) 


0.66(1) 


0.78(4) 


0.69(2) 


0.50(1) 


0.78(3) 




Kcr( K cr) 


0.78(1) 


0.66(1) 


0.78(1) 


0.67(1) 


0.47(1) 


0.78(2) 




N-BPT 


0.74 


0.67 


0.75 


0.73 


0.58 


0.73 




TI-BPT 


0.69 


0.62 


0.71 


0.69 


0.51 


0.68 


6.2 (N f = 0) 




0.82(1) 


0.71(1) 


0.80(1) 


0.71(1) 


0.54(1) 


0.79(3) 




N-BPT 


0.78 


0.73 


0.79 


0.77 


0.65 


0.77 




TI-BPT 


0.72 


0.64 


0.73 


0.71 


0.55 


0.71 



Table 3: Values of the renormalisation constants in the RI-MOM scheme at the scale /i = 1/a for the 
two partially quenched simulations at (i — 5.8 (a -1 = 3.2(1) GeV) and 5.6 (a -1 = 2.4(2) GeV) and for 
the quenched simulation at (3 = 6.2 (a -1 = 3.0(1) GeV). The quoted errors are statistical only. The 
non-perturbatively determined RCs are confronted to their perturbative value (BPT), where "TJ" and "JV" 
stand for tie two types of boosting the bare lattice coupling as discussed in the text. Definition of Z q is the 
same one as in eq. (5) of rcf. [13]. 



The corresponding results for the RCs are then referred to as N-BPT and TI-BPT respectively. 
From the results in table 3 we see that the perturbative estimates of Z$ t p are larger than the non- 
perturbative ones. In other words the non-perturbatively computed Z^ 1 = Za/Zp (Z^ WI = 
1/Zs), at (3 = 5.8, is larger by about 15% (6%) than its BPT counterpart. These differences 
directly propagate to the lattice determination of the renormalised quark masses. As we shall see 
in the following this effect explains the differences between our non-perturbatively renormalised 
results and those obtained by the CP-PACS and JLQCD Collaborations [2,3], where the TI-BPT 
has been used. Furthermore, our results agree with those obtained in refs. [8,9], in which the 
non-perturbative renormalisation has been implemented. 



5 Physical light quark masses 

In order to determine the physical values of the light quark masses, we need to investigate the 
dependence of the pseudoscalar meson masses squared on the valence and sea quark masses. From 
the left plot in fig. 4 we see that: (i) the dependence on the valence quark mass is linear, (ii) the 
slopes do not depend on k s , and (iii) the intercepts do depend on k s and they are non-zero even 
when the sea quark mass is sent to zero. We can therefore fit to the form 

M 2 p (k s , k v1 ,k V2 ) = Co(k s ) + Ci m$ WI [\ + C 2 mf WI + C 3 m£ WI ] , (21) 
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y 



0.07 - 
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0.02 0.024 0.028 0.032 
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NN 
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K =0.1541 



0.030 



0.035 



Figure 4: Dependence of the pseudo-scalar meson masses on the valence quark masses m^ WI (left) and of 
the ratios (Mp — Co)/m^ WI (cf. eq. (21)) on the sea quark masses mf WI . The results refer to the case of 
the partially quenched simulations at j3 = 5.8. 



.AWI 



where the valence and sea AWI quark masses are defined in eq. (10). In this notation vn/ v 
m AWI (k s , k Vi , k V2 ), and m AWI = m wi (k s , k s , k s ). Obviously, the fit forms (16) and (21) differ 
not only in the use of differently defined quark masses (AWI vs. VWI), but also in the presence 
of the constant term Co(n s ) in eq. (21). This term accounts for the 0(a)-discretisation effects 
present in either the pseudoscalar meson masses or the AWI quark masses or in both. Clearly, 
discretisation effects in the pseudoscalar meson masses also affect the fit to eq. (16), i.e., the 
determination of the critical hopping parameters K cr (K s ). 

As before, we find that the values of the coefficients of the quadratic terms, C2 and C3, in 
eq. (21) are consistent with zero [C2 = 0.6(18), and C3 = 0.08(63)], which can be also seen from 
the right plot in fig. 4. Therefore we may set C2 = C3 = 0, and from the fit to eq. (21) we obtain 
Ci = 2.05(3), and 

C (k s ) = {0.0007(18)o.i535 , 0.0023(20)o.i538 , 0.00 46(15) .i 5 40, -0.0039(25) .i 5 4i} • (22) 

Our partially quenched estimates of the average up/down (m u d) and of the strange (m s ) quark 
masses are made by substituting on the l.h.s. of eqs. (16) and (21) the physical pion and kaon 
masses and using the lattice spacing value given in either eq. (13) or eq. (15). In other words, to 
get the quark masses in physical units, we solve 

(ml) phys - 



(m 2 K ) phys - 



h x a 1 x m ud , 



Qi -1 

— x a 
2 



x (m s + m ud ) 



(23) 



in the VWI case, and the analogous expressions, with the coefficient Q\ substituted by C\, in 
the AWI case. Finally the bare quark masses are multiplied by the corresponding RI-MOM mass 
renormalisation constants computed at hq = 1/a, already listed in table 3. 

The conversion to the reference MS scheme, at \x = 2 GeV, is made by using the perturbative 
expressions known up to N 3 LO accuracy [22,27], and A^" 2) = 0.245(16)(20) MeV [28]. We 
finally obtain 

mf (2 GeV) = 120(7) MeV , 



VWI 
AWI 



m ud 
m ud 



(2 GeV) = 5.1(4) MeV 
(2 GeV) = 4.3(4) MeV 



rn 



MS, 



2 GeV) = 101(8) MeV. 



(24) 
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Notice that our VWI result is in very good agreement with the value m^ s (2 GeV) = 119(5)(8) MeV 
obtained by the QCDSF-UKQCD collaboration who uses the VWI quark mass definition and the 
RI-MOM non-perturbative renormalisation [8]. Notice also that our AWI result agrees with the 
value obtained by the Alpha collaboration, mf s (2 GeV) = 97(22) MeV [9]. 

5.1 Systematic uncertainties 

o Discretisation I : From quenched studies we learned that a discrepancy between the quark 
mass values obtained by using two a priori equivalent definitions, namely the VWI and 
AWI definitions, is due to discretisation errors. Such a difference has also been observed in 
other lattice studies with Nf = 2 [2,3]. In particular, in ref. [2] this difference is shown to 
diminish as the lattice spacing is reduced so that, in the continuum limit, the quark masses 
defined in the two ways converge to the same (unique) value. From that study it became 
clear that discretisation errors are significantly larger in the VWI determination of the quark 
masses which is why the JLQCD collaboration choose to quote their AWI masses as final 
results, while the difference between the VWI and AWI results is included in the asymmetric 
systematic error. We follow the same procedure here and we obtain 

mjf (2 GeV) = 4.3(4)(±°- 8 ) MeV , mf(2 GeV) = 101(8)(^ 9 ) MeV . (25) 

o Discretisation II : The results reported in eq. (24) refer to the lattice with coupling (3 = 5.8 
(and volume 24 3 x 48). From our data produced at (3 = 5.6 (see table 1) we obtain 

mjf (2 GeV) = 4.2(3)(+°' 7 ) MeV , mf(2 GeV) = 101(6)(i^ 8 ) MeV , (26) 

thus indistinguishable from those obtained at (5 = 5.8. This makes us more confident that 
the discretisation error attributed to our results in eq. (25) is realistic. 

o Lattice spacing: To obtain the results quoted in eq. (25) we used the lattice spacing (13) 
fixed by m p ^ s ' . These results would become larger by 6% if we used a" 1 (15). We will add 
this difference to our final error budget. 

o Renormalisation constants : As discussed in the text the determination of the mass renor- 
malisation factors is sensitive to discretisation errors. They are especially important when 
the unimproved Wilson quark action is used, which is the case with our study. To check for 
these effects we determined Zp/Zs by using the hadronic Ward identities (hWI) [29] (see 
eqs. (16,17) of ref. [13]), obtaining 

(Zp/Z s ) hW1 = 0.67(1) . (27) 

Compared to the results in table 3, where at (3 = 5.8 we find (Zp /Zs) Rl ~ MOM = 0.75(1), this 
result is about 11% smaller. A similar effect is seen at (3 = 5.6. We observe that, either one 
attributes this effect to Zp or to Zs, {Zp/Zs) hWl always goes in the direction of decreasing 
the AWI- VWI mass difference. Therefore, this 11% uncertainty, being already enclosed in the 
19% systematic error due to the AWI- VWI mass difference, will not be added. We finally 
notice that the result obtained for Zp/Zs by using the numerical stochastic perturbation 
theory to 4-loops, i.e., (Zp/Zs) NSPT = 0.77(1) [30], is more consistent with our larger NPR- 
ed (i.e., RI-MOM) value. 



11 



o Finite volume : Our data at (3 = 5.8 and 16 3 x 48 (see appendix 1) indicate large finite volume 
effects, when compared to our main data-set with (3 = 5.8 and 24 3 x 48. As an illustration 
we show in fig. 5 the effective mass plots obtained in the two volumes at the same value of 
the sea quark mass. Such a sensitivity to the finiteness of the lattice box has already been 
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V= 
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t 

Figure 5: Pseudoscalar and vector effective mass plots at two different volumes at (3 — 5.8 and with 
k s = 0.1538. Valence quark masses are equal to that of the sea quark. The physical volumes correspond to 
(1.5 fm) 3 and (1 fm) 3 respectively. 

observed in ref. [31], where it was shown that, for the range of sea quark masses used in our 
study, finite volume effects are large at L ~ 1 fm but become insignificant when working on 
lattices with L > 1.5 fm. We will rely on that conclusion and will assume that finite volume 
effects are negligible with respect to our statistical and other sources of systematic errors. 
For a more refined estimate of these effects calculations on lattices with larger volumes are 
required. 

We add the above systematic errors algebraically and, as a final result, we quote 



MS 

ud 



[2 GeV) = 4.3 ± 0.4+*' 1 MeV 



mf s {2 GeV) = 101 ± 8±f MeV 



N f = 2, (3 = 5.8 
RI-MOM 



(28) 



which are the values already given in the abstract and the introduction of the present paper. 

5.2 Impact of non-perturbative renormalisation and quenching effects 

Before concluding this section we should compare our main results for the light quark masses 
quoted in eq. (28) with those obtained by evaluating the quark mass RCs using one-loop tadpole 
improved BPT: 



m 



m 



MS, 
ud 1 

MS, 



2 GeV) = 3.7±0.3l^' 3 MeV 



2 GeV) 



±7^° MeV 



Nf = 2, [3 = 5.8 
TI-BPT 



(29) 



These results are in good agreement with those available in the literature obtained from lattice stu- 
dies with Nf = 2 in which renormalisation is implemented perturbatively, namely m^ s (2GeV) = 
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88(+ 4 )MeV [2], and m^ s (2GeV) = 84.5(+} 2 7 °) MeV [3]. Therefore, with Wilson-like dynamical 
fermions, the non-perturbative renormalisation increases the values of the quark masses. When 
comparing our results with those of refs. [2, 3], it should be also noted that the latter are affected 
by discretization errors of different size, having been produced by using 0(a)-improved actions 
but at larger values of lattice spacings. An additional difference, to be kept in mind in such a 
comparison, is that the results of ref. [2] have been obtained after performing an extrapolation to 
the continuum limit. 

Finally, when compared with our results obtained in the quenched simulation at (3 = 6.2, 



on a lattice similar in size and resolution and by using the same (unimproved) Wilson quark action, 
we see no evidence for any effect that could be attributed to the presence of the dynamical quarks. 
For a clearer conclusion concerning this point one should work with lighter sea quark masses which 
is one of the main goals of our future simulations. 

6 Conclusions 

In this paper we have presented our results for the light quark masses obtained from a numerical 
simulation of QCD on the lattice with Nf = 2 degenerate dynamical quarks, with masses covering 
the range 3/4 < m q /m p s hys < 3/2 (with respect to the physical strange quark mass). Our main 
values for the strange and for the average up/down quark masses are obtained on a lattice with 
spacing a ~ 0.06 fm and spatial volume (1.5 fm) 3 , and by using the axial Ward identity. An 
important feature of our calculation is that we renormalise the quark masses non-perturbatively 
in the RI-MOM scheme, which is then converted to the MS scheme by using the known 4-loop 
perturbative formulae. We show that 

- within our statistical accuracy and with relatively heavy dynamical quarks, our unquenched 
results are fully consistent with the quenched ones; 4 

- the non-perturbative renormalisation leads to resulting quark masses larger than those renor- 
malised by using 1-loop perturbation theory (even if tadpole improved); 

- the systematic errors are dominated by lattice discretisation artifacts, which can be cured by 
implementing 0(a)-improvement and/or working with several small lattice spacings followed 
by an extrapolation to the continuum limit. 

The finiteness of the lattice volume is expected not to be an important source of systematic errors 
for the sea quark masses explored in our simulation. Further decrease of the dynamical quark 
masses would necessitate, however, to work with larger physical lattice volumes. Our final results 
are given in eq. (28). 
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Appendix 1: Meson masses from two simulations with larger cou- 
pling (f3 = 5.6) and smaller volume (16 3 x 48) 



In this appendix we provide the reader with two sets of results. We first give our values equivalent 
to the ones presented in table 2 but for (3 = 5.6. Then, for the reader to better appreciate the size 
of finite volume effects at j3 = 5.8, we also tabulate the pseudoscalar and vector meson masses on 
the spatial volumes 16 3 and 24 3 in the fully unquenched situations, i.e., with 





K V1 


K V2 


M P 


My 




0.1560 


0.1560 


0.1560 


0.441(5) 


0.541(6) 


0.0669(4) 




0.1560 


0.1575 


0.398(6) 


0.513(8) 


0.0548(5) 




0.1560 


0.1580 


0.383(7) 


0.503(8) 


0.0509(5) 




0.1575 


0.1575 


0.351(7) 


0.482(9) 


0.0430(5) 




0.1575 


0.1580 


0.334(8) 


0.471(10) 


0.0392(6) 




0.1580 


0.1580 


0.317(9) 


0.459(10) 


0.0355(6) 


0.1575 


0.1560 


0.1560 


0.379(3) 


0.456(8) 


0.0515(4) 




0.1560 


0.1575 


0.333(4) 


0.417(13) 


0.0394(4) 




0.1560 


0.1580 


0.316(4) 


0.403(16) 


0.0355(4) 




0.1575 


0.1575 


0.280(4) 


0.374(21) 


0.0276(3) 




0.1575 


0.1580 


0.261(4) 


0.357(27) 


0.0238(3) 




0.1580 


0.1580 


0.241(4) 


0.340(35) 


0.0200(3) 


0.1580 


0.1560 


0.1560 


0.368(3) 


0.468(11) 


0.0466(4) 




0.1560 


0.1575 


0.324(5) 


0.443(15) 


0.0348(4) 




0.1560 


0.1580 


0.308(6) 


0.436(18) 


0.0311(4) 




0.1575 


0.1575 


0.271(7) 


0.419(21) 


0.0234(5) 




0.1575 


0.1580 


0.250(9) 


0.419(27) 


0.0199(5) 




0.1580 


0.1580 


0.227(11) 


0.424(37) 


0.0164(5) 



Table 4: Vaiues of pseudoscalar and vector meson masses, as well as of the AWI quark masses m^ L (a) = 

\ \m v i{a) + m V 2(a)] AWl obtained by using eq. (10), as obtained from the simulation at /3 = 5.6. All results 
are given in lattice units. 





M p 24) M$*> Ap 


M^ M^ Ay 


0.1535 
0.1538 
0.1540 
0.1541 


0.262(4) 0.299(6) 0.14(4) 
0.236(4) 0.272(10) 0.15(4) 
0.221(3) 0.257(8) 0.17(4) 
0.182(7) 0.251(9) 0.38(4) 


0.348(8) 0.402(13) 0.16(5) 
0.321(5) 0.356(17) 0.11(6) 
0.321(9) 0.347(23) 0.08(7) 
0.291(13) 0.328(31) 0.13(8) 



Table 5: Pseudoscalar and vector meson masses, in lattice units, calculated at (3 = 5.8 on the lattices 
24 3 x 48 and 16 3 x 48, respectively. Finite volume effects are quantified by A Py = [M { py - M^)/M^ . 
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Appendix 2: The scale parameter Rq = ro/a 



To determine the Sommer scale parameter Rq [20] we computed the static quark potential V(R) 
which is extracted from the time dependence of the rectangular Wilson loops W(R,t), where R is 
the spatial separation between the static charges, i.e., 



W(R,t) = C(R)e-nm ^ V(R] = &V- (B,*) S lim log ( jjggjL ) 



(31) 



Large time separations are reached after using the so-called HYP (hypercubic blocking) procedure 
proposed in ref. [33] which helps in improving the signal for the Wilson loops. In this way we find 
a stable signal for V(R) for ig t € [8,13]. At intermediate spatial separations between the static 
charges, we then fit the data (shown in fig. 6) to the form 



V(R) = V + aR--- g5V(R) , 



(32) 



where we also account for SV(R), the perturbatively estimated discretisation correction to the 
Coulomb potential [34]. The Sommer scale Rq is defined as 



R 



2 dV(R)\ 



dR 



) 



1.65. 



(33) 



R=Ro 



The fitting window consistent with the form (32) is found for R £ [2, 7] after which the statistical 



0.5 



0.4 



0.3 - 



0.2 - 



0.1 



Ok =0.1535 

sea 

□ k =0.1538 

sea 

K =0.1540 

sea 

<K =0.1541 



R 



10 



Figure 6: Static quark potential obtained from the simulation with Nf — 2 at (3 — 5.8 on the lattice 
24 3 x 48. 



quality of the data for V{R) rapidly deteriorates. From the fit of our data to the form (32) in that 
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window, we obtain 



= 5.8: k s = {0.1535,0.1538,0.1540,0.1541}, 

R = {7.52(17), 7.70(9), 7.78(19), 8.10(20)}, 

yft = {0.155(4), 0.151(2), 0.150(4), 0.144(4)}, (34) 

e = {0.289(10), 0.302(5), 0.294(7), 0.288(13)}, 

V = {0.260(7), 0.268(3), 0.264(5), 0.262(9)}, 

g = {0.033(14), 0.036(9), 0.035(7), 0.015(15)}, 

ordered as the values of n s . Our results for Rq obtained at the same value of (3 but on the smaller 
volume (16 3 x 48) agree with those written in eq. (34), within the statistical errors. 

At (3 = 5.6, we find stability for tg t 6 [3, 5], and then from the fit to eq. (32) in R £ [2, 7], we 
have 



(3 = 5.6: k s = {0.1560,0.1575,0.1580}, 

R = {5.15(5), 5.72(10), 5.98(8)}, 

v 7 ^ = {0.225(3), 0.202(4), 0.195(3)}, (35) 

e = {0.305(11), 0.317(8), 0.308(6)}, 

V = {0.262(9), 0.279(7), 0.275(5)}, 

g = {0.009(11), 0.022(7), 0.020(4)}, 

in good agreement with ref. [31]. 
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